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Abstract — We propose a new fast algorithm for solving one 
of the standard approaches to ill-posed linear inverse problems 
(IPLIP), where a (possibly non-smooth) regularizer is minimized 
under the constraint that the solution explains the observations 
sufficiently well. Although the regularizer and constraint are 
usually convex, several particular features of these problems 
(huge dimensionality, non-smoothness) preclude the use of off- 
the-shelf optimization tools and have stimulated a considerable 
amount of research. In this paper, we propose a new efficient 
algorithm to handle one class of constrained problems (often 
known as basis pursuit denoising) tailored to image recovery 
applications. The proposed algorithm, which belongs to the family 
of augmented Lagrangian methods, can be used to deal with a 
variety of imaging IPLIP, including deconvolution and recon- 
struction from compressive observations (such as MRI), using 
either total-variation or wavelet-based (or, more generally, frame- 
based) regularization. The proposed algorithm is an instance of 
the so-called alternating direction method of multipliers, for which 
convergence sufficient conditions are known; we show that these 
conditions are satisfied by the proposed algorithm. Experiments 
on a set of image restoration and reconstruction benchmark 
problems show that the proposed algorithm is a strong contender 
for the state-of-the-art. 



I. Introduction 

A. Problem Formulation 

Image restoration/reconstruction is one of the earliest and 
most classical linear inverse problems in imaging, dating back 
to the 1960's [3]. In this class of problems, a noisy indirect 
observation y, of an original image x, is modeled as 

y = Bx + n, 

where B is the matrix representation of the direct operator and 
n is noise. As is common, we adopt the vector notation for 
images, where the pixels on an M x N image are stacked into 
X an (iVAf)-vector in, e.g., lexicographic order In the sequel, 
we denote by n = MN the number of elements of x, thus 
X G M", while y G M™ {m and n may be different). 

In the particular case of image deconvolution, B is the 
matrix representation of a convolution operator This type of 
model describes well several physical mechanisms, such as 
relative motion between the camera and the subject (motion 
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blur), bad focusing (defocusing blur), or a number of other 
mechanisms [7]. 

In more general image reconstruction problems, B repre- 
sents some linear direct operator, such as tomographic projec- 
tions (Radon transform), a partially observed (e.g., Fourier) 
transform, or the loss of part of the image pixels. 

The problem of estimating x from y is called a linear 
inverse problem (LIP); for most scenarios of practical interest, 
this is an ill-posed LIP (IPLIP), i.e., matrix B is singular 
and/or very ill-conditioned. Consequently, this IPLIP requires 
some sort of regularization (or prior information, in Bayesian 
inference terms). One way to regularize the problem of es- 
timating X, given y, consists in a constrained optimization 
problem of the form 



(x) subject to ||Bx — y|l2 < e, 



(1) 



where : R" M = M U {— oo, +00} is the regularizer or 
regularization function, and e > a parameter which depends 
on the noise variance. In the case where 0(x) = ||x||i, the 
above problem is usually known as basis pursuit denoising 
(BPD) [16]. The so-called basis pursuit (BP) problem is the 
particular case of (1) for e = 0. 

In recent years, an explosion of interest in problems of the 
form (1) was sparked by the emergence of compressive sensing 
(CS) [13], [22]. The theory of CS provides conditions (on 
matrix B and the degree of sparseness of the original x) under 
which a solution of (1), for (/)(x) — ||x||i, is an optimal (in 
some sense) approximation to the "true" x. 

In most signal/image recovery and CS problems, the best 
results are obtained with non-smooth regularizers; typical 
examples are the total variation (TV) [13], [47] and £1 
(0(x) = ||x||i) norms. 

B. Analysis and Synthesis Formulations 

In a frame-based representation, the unknown image x 
can be represented as a linear combination of the elements 
of some frame, i.e., x = W/3, where /3 S W^, and the 
columns of the nxd matrix W are the elements of a wavelet' 
frame (an orthogonal basis or an overcomplete dictionary). The 
coefficients of this representation are then estimated from the 
noisy image, under one of the well-known sparsity inducing 
regularizers, such as the ti norm (see [8], [21], [24], [26], 

'We will use the generic term "wavelet" to mean any wavelet-like multi- 
scale representation, such as "curvelets", "beamlets", "ridgelets". 
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[29], and further references therein). Formally, this leads to 
the following constrained optimization problem: 

3 = argimn^(/3) subject to ||BW/3 - y||2 < e. (2) 

This formulation is referred to as the synthesis approach [25], 
[49], since it is based on a synthesis equation: x is synthesized 
from its representation coefficients (x = W/3) which are the 
object of the estimation criterion. Naturally, the estimate of 
X is X = W3. Of course, (2) has the form (1) with BW 
replacing B. 

An alternative formulation applies a regularizer directly 
to the unknown image, leading to criteria of the form (1), 
usually called analysis approaches, since they are based on 
a regularizer that analyzes the image itself, rather than the 
coefficients of a representation thereof. Arguably, the best 
known and most often used regularizer in analysis approaches 
to image restoration is the total variation (TV) norm [15], [47]. 

Wavelet-based analysis approaches are also possible and 
have the form 

min(/)(Px) subject to ||Bx — y||2<£, (3) 

X 

where P is some hnear operator (a matrix) [25]. In this paper, 
we always assume that P is the analysis operator associated 
with 1-tight (Parseval) frame, thus P^P = I [41]. 

C. Previous Algorithms 

If the regularizers are convex, problems (l)-(3) are convex, 
but the very high dimension (at least > 10'', often S> 10^) 
of x and y precludes the direct application of off-the-shelf 
optimization algorithms. This difficulty is further amplified by 
the following fact: for any problem of non-trivial dimension, 
matrices B, W, or P cannot be stored expUcitly, and it is 
costly, even impractical, to access portions (lines, columns, 
blocks) of them. On the other hand, matrix-vector products 
involving these matrices (or their conjugate transposes, denote 
by i')^) can be computed quite efficiently. For example, if the 
columns of W contain a wavelet basis or a tight wavelet frame, 
any multiplication of the form Wv or W^v can be performed 
by a fast transform algorithm [41]. Similarly, if B represents 
a convolution, products by B or B^ can be performed with 
the help of the fast Fourier transform (FFT). These facts have 
stimulated the development of special purpose methods, in 
which the only operations involving matrices are matrix- vector 
products. 

Most state-of-the-art methods for dealing with linear inverse 
problems, under convex, non-smooth regularizers (namely, TV 
and £i), consider, rather than (1), the unconstrained problem 

mini||Bx-y||i+T.^(x), (4) 

where r G R+ is the so-called regularization parameter. Of 
course, problems (1) and (4) are equivalent, in the following 
sense: for any s > such that problem (1) is feasible, a 
solution of (1) is either the null vector, or else is a solution 
of (4), for some r > [30], [46]. For solving problems of 
the form (4), some of the state-of-the-art algorithms belong to 
the iterative shrinkage/thresholding (1ST) family [18], [21], 



[29], [35], and its two-step (TwIST [9] and FISTA [5]) 
and accelerated (SpaRSA [53]) variants. These methods were 
shown to be considerably faster than earlier methods, including 
ll_ls [37] and the codes in the ii-magic-^ and the SparseLab^ 
toolboxes. 

A key ingredient of most of these algorithms is the so- 
called shrinkage/thresholding/denoising function, which is the 
Moreau proximal mapping of the regularizer [18]. Formally, 
this function : ^ M" is defined as 

*r^(y) = argmm i||x - y\\l + t(/)(x). (5) 

Notice that if ^ is proper and convex, the function being 
minimized is proper and strictly convex, thus the minimizer 
exists and is unique making the function well defined [18]. For 
some choices of (j), the corresponding ^r<j) have well known 
closed forms. For example, if (/)(x) = ||x||i, then *T-^(y) = 
soft(y, r), where soft(-, r) denotes the component-wise appli- 
cation of the soft-threshold function y i— > sign(?/) max{|?/| — 
r,0}. 

In [27], [2], we proposed a new algorithm called split aug- 
mented Lagrangian shrinkage algorithm (SALSA), to solve 
unconstrained optimization problems of the form (4) based 
on variable splitting [19], [52]. The idea is to transform the 
unconstrained problem (4) into a constrained one via a variable 
splitting "trick", and then attack this constrained problem 
using an augmented Lagrangian (AL) method [44]. AL is 
known to be equivalent to the Bregman iterations recently 
proposed to handle imaging inverse problems (see [54] and 
references therein). We prefer the AL perspective, rather than 
the Bregman iterative view, as it is a more standard and 
elementary tool (covered in most optimization textbooks). On 
several benchmark experiments (namely image deconvolution, 
recovery of missing pixels, and reconstruction from partial 
Fourier observations) using either frame-based or TV-based 
regularization, SALSA was found to be faster than the previous 
state-of-the-art methods FTSTA [5], TwIST [9], and SpaRSA 
[53]. 

Although it is usually easier to solve an unconstrained prob- 
lem than a constrained one, formulation (1) has an important 
advantage: parameter e has a clear meaning (it is proportional 
to the noise standard deviation) and is much easier to set than 
parameter r in (4). Of course, one may solve (1) by solving 
(4) and searching for the "correct" value of r that makes (4) 
equivalent to (1). Clearly, this is not efficient, as it involves 
solving many instances of (4). Obtaining fast algorithms for 
solving (1) is thus an important research front. 

There are few efficient algorithms to solve (l)-(3) in an 
image recovery context: x and y of dimension > 10^ (often 
> 10^), B representing an operator, and (j) a convex, non- 
smooth function. A notable exception is the recent SPGLl 
[51], which (as its name implies) is specifically designed for 
£i regularization. As shown in [51], the methods for solving 
(1) available in the £i-magic package are quite inefficient for 
large problems. General purpose methods, such as the SeDuMi 

^Available at http: //www. 11-magic . org 
'Available at http: // spar selab. Stanford, edu 
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package"*, are simply not applicable when B is not an actual 
matrix, but an operator. 

Another efficient algorithm for solving problems of the 
form (1) is the recently proposed NESTA [6], which is 
based on Nesterov's first-order methods which achieve an 
optimal convergence rate by coupling smoothing techniques 
with an improved gradient method [42], [43]. NESTA allows 
for minimizing either the £i or TV norm, and also allows 
using synthesis or analysis formulations. Nesterov's first-order 
method was also recently adopted in [20] to perform TV- 
regularized image denoising, deblurring, and inpainting. 

The Bregman iterative algorithm (BIA) was recently pro- 
posed to solve (1) with £ = 0, but is not directly applicable 
when e > [54]. To deal with the case of £ > 0, it was 
suggested that the BIA for £ = is used and stopped when 
||Bx-y||2 < £ [11], [54]. Clearly, that approach is not 
guaranteed to find a good solution, since it depends strongly 
on the initialization; e.g., if the algorithm starts at a feasible 
point, it will immediately stop, although the point may be far 
from a minimizer of (j). 

D. Proposed Approach 

In this paper, we introduce an algorithm for solving opti- 
mization problems of the form (1). The original constrained 
problem (1) is transformed into an unconstrained one by 
adding the indicator function of the feasible set, the elhpsoid 
{x : ||Bx — y|| < s}, to the objective in (1). The resulting 
unconstrained problem is then transformed into a different 
constrained problem, by the application of a variable split- 
ting operation; finally, the obtained constrained problem is 
dealt with using the alternating direction method of multipli- 
ers (ADMM) [23], [31], [32], which belongs to the family 
of augmented Lagrangian (AL) techniques [44]. Since (as 
SALSA), the proposed method uses variable sphtting and AL 
optimization, we call it C-SALSA (for constrained-SAl^SA). 

The resulting algorithm is more general than SPGLl, in 
the sense that it can be used with any convex regularizer 4> 
for which the corresponding Moreau proximity operator (see 
[18]), has closed form or can be efficiently computed. In this 
paper, we will show examples of C-SALSA where x is an 
image, (j> is the TV norm [47], and ^t4> is computed using 
Chambolle's algorithm [14]. Another classical choice which 
we will demonstrate is the £i norm, which leads to (y) = 
soft(y,r). 

C-SALSA is experimentally shown to efficiently solve im- 
age recovery problems, such as MRI reconstruction from CS- 
type partial Fourier observations using TV regularization, and 
image deblurring using wavelet-based or TV regularization, 
faster than SPGLl and NESTA. 

E. Organization of the Paper 

The paper is organized as follows. Section 11 describes the 
basic ingredients of C-SALSA: variable splitting, augmented 
Lagrangians, and the ADMM. Section III contains the deriva- 
tion leading to C-SALSA. Section IV reports experimental 

''Available at http://sedumi.ie.lehigh.edu 



results, and Section V ends the paper with a few remarks and 
pointers to future work. 

II. Basic Ingredients 
A. Variable Splitting 

Consider an unconstrained optimization problem 



min/i(u) + /2(Gu), 



(6) 



where G e M'*''", /i : M" ^ M, and /s : M*^ ^ M. Variable 
splitting (VS) is a simple procedure that consists in creating 
a new variable, say v, to serve as the argument of /2, under 
the constraints that v = Gu, i.e.. 



mm 



/i(u)-F/2(v), subject to v = Gu. (7) 



The rationale behind VS is that it may be easier to solve 
the constrained problem (7) than it is to solve its equivalent 
unconstrained counterpart (6). 

VS has been recently used in several image processing 
applications. A VS method was used in [52] to obtain a 
fast algorithm for TV-based restoration. Variable splitting was 
also used in [10] to handle problems involving compound 
regularizers. In [10] and [52], the constrained problem (7) is 
attacked using a quadratic penalty approach, i.e., by solving 



Q. 

min /i(u) + /2(v) + - ||Gu • 

ueM^jveK"* 2 



112) 



(8) 



by alternating minimization with respect to u and v, while 
slowly taking a to very large values (a continuation process), 
to force the solution of (8) to approach that of (7), which in 
turn is equivalent to (6). The rationale of these methods is that 
each step of this alternating minimization may be much easier 
than the original unconstrained problem (6). The drawback 
is that as a grows, the intermediate minimization problems 
become increasingly ill-conditioned, thus causing numerical 
problems [44]. 

A similar VS approach underlies the recently proposed split- 
Bregman methods [33]; however, instead of using a quadratic 
penalty technique, those methods attack the constrained prob- 
lem directly using a Bregman iterative algorithm [54], which is 
known to be equivalent to the augmented Lagrangian method 
[34], [50], [54]. 

B. Augmented Lagrangian 

Consider the constrained optimization problem with linear 
equahty constraints 

min Eiz) subject to Az - b = 0, (9) 

zeM" 

where b G Rp and A G M^^", i.e., there are p linear 
equahty constraints. The so-called augmented Lagrangian for 
this problem is defined as 



£a(z, a, = E{z,) + X^{h - Az) 



Az-b||2, (10) 



where A € Rp is a vector of Lagrange multipliers and /x > is 
called the AL penalty parameter [44]. The so-called augmented 
Lagrangian method (ALM) [44], also known as the method 
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of multipliers (MM) [36], [45], iterates between minimizing 
£^(z, A, /Lt) with respect to z, keeping A fixed, and updating 
A, until some convergence criterion is satisfied: 

Algorithm ALM/MM 

1. Set fc = 0, choose /z > and Aq. 

2. repeat 

3. Zfe+i G argminj./:A(z, Aft,/i) 

4. Afe+i = Afe + (U(Azfe+i - b) 

5. + 1 

6. until some stopping criterion is satisfied. 

It is also possible (and even recommended) to update the 
value of jj, in each iteration [4], [44]. However, unlike in the 
quadratic penalty approach, the ALM/MM does not require 
to be taken to infinity to guarantee convergence to the solution 
of the constrained problem (9). 

After a straightforward complete-the-squares procedure, 
the terms added to -E(z) in the augmented Lagrangian 
Afe, fx) can be written as a single quadratic term (plus a 
constant independent of z, thus irrelevant to the ALM/MM), 
leading to the following alternative form of the algorithm 
(which makes clear its equivalence with the Bregman iterative 
method [54]): 

Algorithm ALM/MM (version II) 



1. 
2. 
3. 
4. 
5. 
6. 



Set k = 0, choose /x > and do. 
repeat 

Zfe+i e argmiiiz E{z) + f ||Az 
dfc+i = dfc - (Azfe+i - b) 

witil some stopping criterion is satisfied. 



k\\2 



C. ALM/MM for Variable Splitting and ADMM 

The constrained problem (7) can be written as (9) by 
defining £(z) = /i(u) + /2(v) and setting 

z= [u'^v'^]'^, b = 0, A=[G -I]. (11) 

With these definitions in place. Steps 3 and 4 of the ALM/MM 
(version II) become 

(ufe+i, Vfe+i) e argmin/i(u) + /2(v) + ^||Gu - v - dk||2 

u,v Z 

dfe+1 = dfe - (Gufc+i - Vfe+i). (12) 

The minimization problem yielding (ufe+i, Vfe+i) is not trivial 
since, in general, it involves a non-separable quadratic term 
and possibly non-smooth terms. A natural approach is to use 
a non-linear block-Gauss-Seidel (NLBGS) technique which 
alternates between minimizing with respect to u and v while 
keeping the other fixed. Of course this raises several questions: 
for a given dfe, how much computational effort should be spent 
in this problem? Does the NLBGS procedure converge? The 
simplest answer to these questions is given in the form of the 
so-called alternating direction method of multipliers (ADMM) 
[23], [31], [32], which is simply an ALM/MM in which only 
one NLBGS step is performed in each outer iteration. 



Algorithm ADMM 

1. Set A: = 0, choose p> Q, vq, do. 

2. repeat 

3. Ufe+i e argminu/i(u) + fUGu- Vfe - dfe||| 

4. Vfe+i G argminv /2(v) + f ||Gufe+i - v - dk\\l 

5. dfe+i = dfe — (Gufc+i — Vfe-|_i) 

6. k^k + l 

1. until some stopping criterion is satisfied. 

For later reference, we now recall the theorem by Eckstein 
and Bertsekas [23, Theorem 8] in which convergence of (a 
generalized version of) ADMM is shown. 

Theorem 1 (Eckstein-Bertsekas, [23]}: Consider problem 
(6), where G has full column rank, and fi and are 
closed, proper, convex functions. Consider arbitrary /z > 
and uo,do,vo G W. Let {rjk > 0, k = 0,1,...} and 
{vk > 0, fc = 0, 1, ...} be two sequences such that 

oo oo 

^^77fe<oo and Vk < oo. 



fc=o 



fe=o 



Consider three sequences {uk G R", k = 0,1,...}, {vfe G 
M'', A; = 0,l,...},and{dfc G R'', k = 0,1, ...} that satisfy 



Ufe+i - argmin/i(u) + ^||Gu- Vfe-dfe||2 
u 2 

Vfe+i - argmin/2(v) + ^||Gufe+i-v-dfc||i 
d/s+i = dfe - (G Ufe+i - Vfe+i). 



< 
< 



Vk 



Then, if (6) has a solution, say u*, then the sequence {ufe} 
converges to u*. If (6) does not have a solution, then at least 
one of the sequences {ufe} or {d^} diverges. 

Notice that the ADMM as defined above (if each step is 
implemented exactly) generates sequences {ufe}, {vfe}, and 
{dfe} that satisfy the conditions in Theorem 1 in a strict sense 
(i.e., with rjk = = 0). The remaining key condition for 
convergence is then that G has full column rank. One of the 
important corollaries of this theorem is that it is not necessary 
to exactly solve the minimizations in lines 3 and 4 of ADMM; 
as long as the sequence of errors are absolutely summable, 
convergence is not compromised. 

The proof of Theorem 1 is based on the equivalence between 
ADMM and the Douglas-Rachford Splitting (DRS) applied to 
the dual of problem (6). The DRS was recently used for image 
recovery problems in [17]. For recent and comprehensive 
reviews of ALM, ADMM, DRS, and their relationship with 
Bregman and split-Bregman methods, see [34], [50]. 

D. A Variant of ADMM 

Consider a generalization of problem (6), where instead of 
two functions, there are J functions, that is. 



min o, (H'-''' u), 



(13) 



where <?, : W' — > M are closed, proper, convex functions, 
and HW G Wi""^ are arbitrary matrices. The mirumization 
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problem (13) can be written as (6) using the following corre- 
spondences: /i = 0, 

H(i) 



G = 



H(-^) 



(14) 



where p = pi 



^pxd 



■ +PJ, and /2 



given by 



(15) 



where v(j) G and v = [(v(i))^, . . . , (v('^))'^]^ G W. 
We now simply apply ADMM (as given in the previous 
subsection), with 





■ 4^^ " 








. 4-" . 


, Vfc = 





Moreover, the fact that /i = turns Step 3 of the algorithm 
into a simple quadratic minimization problem, which has a 
unique solution if G has fuU column rank: 

argmin||Gu-Cfc||2 = (G^G)"'g^C„ 



(16) 



■i=i 



k ' 



i=i 



a) a) a) 

where Ck =^k + dfe (and, naturally, Ck — ^k + '^fc ) 
the second equality results from the particular structure of G 
in (14). 

Furthermore, our particular way of mapping problem (13) 
into problem (6) allows decoupling the minimization in Step 
4 of ADMM into a set of J independent ones. In fact, 

Vfe+i <- argmin/2(v) + ^||Gufe+i - v-dk\\l (17) 

can be written as 



L fc+i J 



+ 



arg „min giiv^'>) + ■ ■ ■ + gjiv^'>)+ 

v' ' .... .v' ' 





■ H(i) ■ 
















Ufe+l - 











Clearly, the minimizations with respect to u^^\ . . . . 
decoupled, thus can be solved separately, leading to 



2 

u(^) 



are 



(i) ■ I \ w 

' arg mm qA'v) H — v 



,(1)||2 



(18) 



for 2 = 1, J, where 

s«=HWu,+i-d«. 

Since this algorithm is exactly an ADMM, and since all the 
functions g^, for j = 1, J, are closed, proper, and convex, 
convergence is guaranteed if G has fuU column rank. Actually, 
this fuU column rank condition is also required for the inverse 



in (16) to exist. Finally, notice that the update equations in 
(18) can be written as 



(19) 



where the 'S'gj/p are, by definition, the Moreau proximal 
mappings of gi/n, gj/fi. 

In summary, the variant of ADMM (herein referred to as 
ADMM-2) that results from the formulation just presented is 
described in the following algorithmic framework. 



Algorithm ADMM-2 
1. 

2. repeat 

3. for i = 1, J 



Set k = 0, choose > 0, v^^^ \ d[^^\ d^-^^ 



do 



Ufe+l 

for i = r, J 

^k+i 



do vli , = (H(')ufe+i 



-40 



9. + 1 

10. until some stopping criterion is satisfied. 



III. Proposed Method 

We now apply the algorithmic framework described in the 

previous section to the basic problem (1) (which includes (2) 
as a special case), as well as the analysis formulation (3). 

A. Problem (1) 

For the constrained optimization problem (1), the feasible 
set is the ellipsoid 



^(e,B,y) = {xe 



|Bx-y||2 < e}, 



(20) 



which is possibly infinite in some directions (since B may 
be singular). Problem (1) can be written as an unconstrained 
problem, with a discontinuous objective. 



where is '■ 
S C M™, 



niin (f>{x) + t£;(£j,y)(Bx), (21) 
— > M denotes the indicator function of set 



i-sis) 



0, ifse5 
+00, if s ^ S. 



(22) 



Notice that E{e.l,y) is simply a closed e-radius EucUdean 
ball centered at y. 

Problem (21) has the form (13) with J = 2 and 



gi = 9 

52 = ''E{e,I,y) ) 

H(i) ^ I 

H(2) = B. 



(23) 
(24) 
(25) 
(26) 



Instantiating ADMM-2 to this particular case requires the 
definition of the Moreau proximal maps associated with gi = 
(j) and 52 = i'E{e,i,y) - Concerning </>, the regularizer, we assume 
that ^T(t>{') (see (5)) can be computed efficiently. This is of 
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course the case of (/)(x) = ||x||i, for which is simply 
a soft threshold. If (p is the TV norm, we may use one 
the fast algorithms available to compute the corresponding 
denoising function [14], [20]. The Moreau proximal map of 
g2 = iE{e,i,y) is defined as 



*^i.(e,.,y)/M(s) = argmin. 



x-s||i, (27) 



which is obviously independent of and is simply the or- 
thogonal projection of s on the closed e-radius ball centered 
at y: 



if l|s -y||2 > e, 
if l|s - y||2 < e. 



(28) 



We are now in a position to instantiate ADMM-2 for solving 
(21) (equivalently (1)). The resulting algorithm, which we call 
C-SALSA-1, is as follows. 



Algorithm C-SALSA-1 

1. SetA; = 

2. repeat 



1. Set fc = 0, choose > 0, v^^\ v^^\ d^^K d^^K 



3. 
4. 

5. 
6. 
7. 
8. 
9. 



U/c+l 



I V 

-1 



(2) 




j(2) 



I + B^B 



.+1 - (ufc+1 - 

1 = *^.(e.I.., 



^(2^ _d(2) 



fe+i 

Ufe+l + V 

Bufe+i + V 



,(2)^ 



(1) 
fe+1 
(2) 
fc+1 



10. until some stopping criterion is satisfied. 

The issue of how to efficiently solve the linear system 
of equations in line 4 of C-SALSA-1 will be addressed in 
Subsection 111-C. 

Convergence of C-SALSA-1 is guaranteed by Theorem 1 
since it is an instance of ADMM with 



G = 



I 
B 



(29) 



which is a full column rank matrix, and both (f) and LE{e,i,y) 
are closed, proper, convex functions. 

Finally, notice that to apply C-SALSA-1 to problem (2) we 
simply have to replace B with BW. 

B. Problem (3) 

Problem (3) can also be written as an unconstrained problem 

min (/>(Px) + tB(e,i,y) (Bx), (30) 

which has the form (13) with J = 2 and 

91 = cj> (31) 

92 = l'E{s,I,y), (32) 

H(i) = P (33) 
H(2) = B. (34) 



The resulting ADMM algorithm, called C-SALSA-2, is 
similar to C-SALSA-1, with only a few minor differences. 

Algorithm C-SALSA-2 

1. Set fc = 0, choose /x > 0, v^^^ v^^K d^o \ d^^\ 

2. repeat 

3. r.=p(v«+d«)+B-(vf)+dl^)) 

4. Ufc+i=( P«P + B«B I rfc 



Ufc + l : 

5. 
6. 

7. dl 



(Pufe+1 - dfe 



(1) 



E(e,I,y) 
(1) 



Bu 



fe+i 



PUfe+l + v^j^li 



1^'^ -d 

fc+1 — k 

o h(2) _h(2) R„ -Ui/2) 

9. + 1 

10. until some stopping criterion is satisfied. 



In this paper, we assume that P is the analysis operator 
of a 1-tight (Parseval) frame, thus P^P = I and Une 4 of 
C-SALSA-2 is similar to line 4 of C-SALSA-1: 



Ufc+i = (l + B^B) ' 



(35) 



The issue of how to efficiently solve this linear system will 
be addressed in the next subsection. 

Since both (j) and iB(e,i,y) are closed, proper, convex func- 
tions, convergence of C-SALSA-2 holds (by Theorem 1) if 



G 



P 
B 



(36) 



is a full column rank matrix. This is of course true if P is 
itself a full column rank matrix, which is the case if P is the 
analysis operator of a tight frame [41]. 



C. Solving (35) 

As mentioned in Subsection 1-C, in most imaging problems 
of interest, it may not be feasible to explicitly form matrix 
B. This might suggest that it is not easy, or even feasible, 
to compute the inverse of (l-|-B^B). However, as shown 
next, in a number of problems of interest, this inverse can be 
computed very efficiently with O(nlogn) cost. 

1) Deconvolution with Analysis Formulation: In the case 
of analysis formulations of the form (1) or (3) to image de- 
convolution problems, matrix B represents a 2D convolution. 
Consequently, matrix B can be factorized as B = U^DU, 
where U is the unitary matrix (U^ = U~^) representing the 
discrete Fourier transform (DFT) and D is diagonal. Thus, 



(B^B + I)-^ =U^ (|D 



(37) 



where |D|^ is the matrix with squared absolute values of the 
entries of D. Since |Dp +1 is diagonal, its inversion cost is 
0{n). Products by U and have O(nlogn) cost, using the 
FFT algorithm. 
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2) Deconvolution with Synthesis Formulation: In this case, 
as seen in Section I-B, we have BW instead of B, and even 
if B is a convolution, BW is not diagonalizable by the DFT. 
To sidestep this difficulty, we assume that W contains a 1- 
tight (Parseval) frame {i.e., W = I). Using the Sherman- 
Morrison- Woodbury (SMW) matrix inversion lemma, 

(W'^B^BW-M)"^= I - W-^B" (bB^ -I- 1) ^BW, 



thus Une 4 of C-SALSA-1 and C-SALSA-2 can be written as 

Ufe+i = (rfc - W^FW Vk) . (39) 

Since B is a convolution, B = U^DU, thus multiplying by F 
corresponds to applying an image filter in the Fourier domain 

F = U-^D* (|D|2 + 1)"^ DU, 

which has O(nlogn) cost, since all the matrices in 
D* (|Dp + 1) D are diagonal and the products by U and 
are carried out via the FFT. The cost of (38) will thus be 
either 0{n log n) or the cost of the products by and W. 

For most tight frames used in image processing, there are 
fast O(nlogn) algorithms to compute the products by 
and W [41]. For example, in the case of translation-invariant 
wavelet transforms, these products can be computed using 
the undecimated wavelet transform with 0(n log n) total cost 
[39 J. Curvelets also constitute a Parseval frame for which 
fast 0(n log n) implementations of the forward and inverse 
transform exist [12]. Yet another example of a redundant 
Parseval frame is the complex wavelet transform, which has 
0(n) computational cost [38], [48]. In conclusion, for a large 
class of choices of W, each iteration of the SALSA algorithm 
has O(nlogn) cost. 

3) Missing Pixels: Image Inpainting: In the analysis prior 
formulation of this problem, the observation matrix B models 
the loss of some image pixels; it is thus an m x n binary 
matrix, with m < n, which can be obtained by taking a subset 
of rows of an identity matrix. Due to its particular structure, 
this matrix satisfies BB^ = I. Using this fact together with 
the SMW formula leads to 

{B"B + iy^ = I- B^ (I + BB^)~^B (40) 

= (41) 

Since B^B is equal to an identity matrix with some zeros 
in the diagonal (corresponding to the positions of the missing 
observations), the matrix in (41) is diagonal with elements 
either equal to 1 or 1/2. Consequently, line 4 of C-SALSA- 
1 and C-SALSA-2 corresponds to multiplying this diagonal 
matrix by rfc, obviously with 0{n) cost. 

In the frame-based synthesis formulation, we have BW 
instead of B. Using the SMW formula yet again, and the facts 
that BB^ = I and WW^ = I, we have 

(W^B-^BW -h I) = I - - W^B^BW. (42) 



As noted in the previous paragraph, A^A is equal to an 
identity matrix with zeros in the diagonal, i.e., a binary 
mask. Thus, the multiplication by W^A^AW corresponds 
to synthesizing the image, multiplying it by this mask, and 
computing the representation coefficients of the result. In 
conclusion, the cost of line 4 of C-SALSA-1 and C-SALSA- 
2 is again that of the products by W and W^, usually 
0{n log n). 

4) Partial Fourier Observations (MRI Reconstruction): 
Finally, we consider the case of partial Fourier observations, 
which is used to model MRI acquisition and has been the focus 
of recent interest due to its connection to compressed sensing 
[13], [40]. In the analysis formulation, B = MU, where M 
is an m X n binary matrix (m < n) again, formed by a subset 
of rows of the identity, and U is the DFT matrix. Due to its 
particular structure, matrix M satisfies MM^ = I; this fact 
together with the matrix inversion lemma leads to 

(B^B + =1-^ U^M^MU, (43) 

where M^M is equal to an identity with some zeros in the 
diagonal. Consequently, the cost of Une 4 of C-SALSA-1 and 
C-SALSA-2 is again that of the products by U and U^, i.e. 
0{n\ogn) using the FFT. 

In the synthesis case, the observation matrix has the form 
MUW. Clearly, the case is again similar to (42), but with 
UW and W^U^ instead of W and W^, respectively. 
Again, the cost of Une 4 of C-SALSA-1 and C-SALSA-2 is 
O(nlogn), if the FFT is used to compute the products by U 
and and fast frame transforms are used for the products 
by W and W^. 

D. Computational Complexity 

As shown in the previous section, the cost of line 4 of 
C-SALSA-1 and C-SALSA-2 is O(nlogn). The other Unes 
of the algorithms simply involve: (a) matrix- vector products 
involving B, W, P, or their conjugate transposes, which have 
0(n log n) cost; (b) vector additions, which have 0{n) cost; 
and (c) the computation of the Moreau proximal maps (lines 
5 and 6 of C-SALSA-1 and C-SALSA-2). In the case of the 
projections on a ball (line 6), it is clear from (28) that the cost 
is 0{n). 

Finally, we consider the computational cost of the Moreau 
proximal map of the regularizer (j) (Une 5 of C-SALSA-1 
and C-SALSA-2). In some cases, this map can be computed 
exactly in closed form; for example, if 0(x) = ||x|| i, then 
is simply a soft threshold and the cost is 0{n). In other cases, 
the Moreau proximal map does not have a closed form solu- 
tion; for example, if (/)(x) = TV(x), the corresponding 
has to be computed using one of several available iterative 
algorithms [14], [20]. Most of these iterative algorithms can 
be implemented with 0{n) cost, although with a factor that 
depends on the number of iterations. In our implementation 
of C-SALSA we use Chambolle's algorithm [14]. 

In summary, for a wide choice of regularizers and frame 
representations, the C-SALSA algorithms have O(nlogn) 
computational complexity. 
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IV. Experiments 

In this section, we report results of experiments aimed at 
comparing the speed of C-SALSA with that of the current state 
of the art methods (that are freely available online): SPGLl^ 
[51], and NESTA*' [6]. 

We consider three standard and often studied imaging 
inverse problems: image deconvolution (using both wavelet 
and TV-based regularization); image restoration from missing 
samples (inpainting); image reconstruction from partial Fourier 
observations, which (as mentioned above) has been the focus 
of much recent interest due to its connection with compressed 
sensing and the fact that it models MRI acquisition [40]. 

All experiments were performed using MATLAB, on a 
Windows XP desktop computer with an Intel Pentium-IV 3.0 
GHz processor and 1.5 GB of RAM. The number of calls 
to the operators B and B^, the number of iterations, CPU 
times, and MSE values presented are the averages values over 
10 runs of each experiment. The number of calls reported for 
each experiment is the average over the 10 instances, with the 
minimum and maximum indicated in the parentheses. Since 
the stopping criteria of the implementations of the available 
algorithms differ, to compare the speed of the algorithms in a 
way that is as independent as possible from these criteria, the 
experimental protocol that we followed was the following: we 
first run one of the algorithms with its stopping criterion, and 
then run C-SALSA until the constraint in (1) is satisfied and 
the MSE of the estimate is below that obtained by the other 
algorithms. 

The value of e in (1) used in all cases was y^m + S^^/ma, 
where m is the number of observations, and a is the noise 
standard deviation. The parameter /i was hand-tuned for fastest 
convergence. 

A. Image Deconvolution with wavelets 

We consider five benchmark deblurring problems [29], sum- 
marized in Table I, all on the well-known Cameraman image. 
The regularizer is <j){f3) = ||/3||i, thus ^^^^ is an element- 
wise soft threshold. We compare C-SALSA against SPGLl 
and NESTA in the synthesis case, and against only NESTA 
in the analysis case, since SPGLl is hardwired with ||x||i as 
the regularizer, and not ||Px||i. Since the restored images are 
visually indistinguishable from those obtained in [29], and the 
SNR improvements are also very similar, we simply compare 
the speed of the algorithms, that is, the number of calls to 
the operators B and B^, the number of iterations, and the 
computation time. 

TABLE I 

Details of the image deconvolution experiments. 



Experiment 


blur kernel 


a2 


1 


9x9 uniform 


0.56^ 


2A 


Gaussian 


2 


2B 


Gaussian 


8 


3A 


h^j = l/{l + P+f) 


2 


3B 




8 



^Available at http : / / www . cs . ubc .ca/labs/scl/spgll 
^Available at http : / /www . acm . caltech . edu/ -nesta 



10* 



lO"* 











SALSA 




■ SPGll 




NESTA 



100 200 300 400 500 600 700 
seconds 

(a) 



10* 





• e 




SALSA 




• SPGll 




NESTA 


I ! „r'"";rr:- 











100 200 300 40U 500 600 700 

seconds 



(b) 

Fig. 1. Image deblurring with wavelets (synthesis prior, redundant Haar 
wavelets), 9x9 uniform blur, cr = 0.56: (a) Evolution of the objective 
function ||x||i over time; (b) quadratic constraint ||AWx — y||2 over time. 



In the first set of experiments, W is a redundant Haar 
wavelet frame with four levels. For the synthesis case, the CPU 
times taken by each of the algorithms are presented in Table 11. 
Table III presents the corresponding results for the case with 
the analysis prior. In the second set of experiments, W is 
an orthogonal Haar wavelet basis; the results are reported in 
Table IV for the synthesis case, and in Table V for the analysis 
case. To visually illustrate the relative speed of the algorithms. 
Figure 1 plots the evolution of the constraint ||Bufe — y||, 
versus time, in experiments 1, for the synthesis prior case, 
with redundant wavelets. 



B. Image Deblurring with Total Variation 

The same five image deconvolution problems listed in Ta- 
ble I were also addressed using total variation (TV) regulariza- 
tion (more specifically, the isotropic discrete total variation, as 
defined in [14]). The corresponding Moreau proximal mapping 
is computed using 5 iterations of Chambolle's algorithm [14]. 

Table VI compares the performance of C-SALSA and 
NESTA, in terms of speed. The evolutions of the objective 
function and the constraint for experiment 1 are plotted in 
Figure 2. 

We can conclude from Tables II, III, IV, V, and VI that, in 
image deconvolution problems, both with wavelet-based and 
TV-based regularization, C-SALSA is almost always clearly 
faster than the fastest of the other competing algorithms. 
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TABLE II 

Image deblurring using wavelets (redundant) - Computational load 



Expt. 


Avg. calls to B,B^ (min/max) 


Iterations 


CPU time (seconds) 




SPGLl 


NESTA 


c-salsa 


SPGLl 


NESTA 


C-SALSA 


SPGLl 


NESTA 


C-SALSA 


1 


1029 (659/1290) 


3520 (3501/3541) 


398 (388/406) 


340 


880 


134 


441.16 


590.79 


100.72 


2A 


511 (279/663) 


4897 (4777/4981) 


451 (442/460) 


160 


1224 


136 


202.67 


798.81 


98.85 


2B 


377 (141/532) 


3397 (3345/3473) 


362 (355/370) 


98 


849 


109 


120.50 


557.02 


81.69 


3A 


675 (378/772) 


2622 (2589/2661) 


172 (166/175) 


235 


656 


58 


266.41 


423.41 


42.56 


3B 


404 (300/475) 


2446 (2401/2485) 


134 (130/136) 


147 


551 


41 


161.17 


354.59 


29.57 



TABLE III 

Image Deconvolution using wavelets (redundant, analysis prior) - Computational load 



Expt. 


Avg. calls to B,B^ (min/max) 


Iterations 


CPU time (seconds) 




NESTA 


C-SALSA 


NESTA 


C-SALSA 


NESTA 


C-SALSA 


1 


2881 (2861/2889) 


413 (404/419) 


720 


138 


353.88 


80.32 


2A 


2451 (2377/2505) 


362 (344/371) 


613 


109 


291.14 


62.65 


2B 


2139 (2065/2197) 


290 (278/299) 


535 


87 


254.94 


50.14 


3A 


2203 (2181/2217) 


137 (134/143) 


551 


42 


261.89 


23.83 


3B 


1967 (1949/1985) 


116 (113/119) 


492 


39 


236.45 


22.38 



TABLE IV 

Image deblurring using wavelets (orthogonal) - Computational load 



Expt. 


Avg. calls to B, B^ (min/max) 


Iterations 


CPU time (seconds) 




SPGLl 


NESTA 


C-SALSA 


SPGLl 


NESTA 


C-SALSA 


SPGLl 


NESTA 


C-SALSA 


1 


730 (382/922) 


13901 (13871/13931) 


494 (424/748) 


298 


3475 


166 


46.64 


622.09 


23.91 


2A 


352 (191/480) 


1322 (1301/1329) 


205 (202/205) 


128 


331 


69 


19.21 


58.29 


10.07 


2B 


207 (128/254) 


1218 (1193/1261) 


123 (115/133) 


87 


305 


42 


12.23 


52.92 


6.35 


3A 


248 (161/320) 


1421 (1413/1433) 


118 (115/121) 


104 


355 


40 


14.98 


58.693 


5.57 


3B 


170 (114/220) 


4408 (4345/4545) 


258 (94/328) 


72 


1102 


87 


9.51 


181.83 


11.93 



TABLE V 

Image deblurring using wavelets (orthogonal, analysis prior) - Computational load 



Expt. 


Avg. calls to B, B^ (min/max) 


Iterations 


CPU time (seconds) 




NESTA 


C-SALSA 


NESTA 


C-SALSA 


NESTA 


C-SALSA 


1 


8471 (8413/8553) 


387 (380/395) 


2118 


117 


300.60 


16.51 


2A 


2463 (2445/2489) 


377 (371/383) 


616 


126 


311.49 


77.75 


2B 


2159 (2097/2253) 


300 (290/317) 


540 


101 


280.35 


59.75 


3A 


2203 (2165/2229) 


153 (149/155) 


551 


52 


282.12 


32.02 


3B 


4710 (4577/4829) 


212 (104/374) 


1178 


59 


167.73 


7.89 



C. MRI Image Reconstruction 

We now consider the problem of reconstructing the 128 x 
128 Shepp-Logan phantom (shown in Figure 3(a)) from a 
limited number of radial lines (22, in our experiments, as 
shown in Figure 3(b)) of its 2D discrete Fourier transform. 
The projections are also corrupted with circular complex 
Gaussian noise, with variance fj^ = 0.5 x 10~^. We use 
TV regularization (as described in Subsection IV-B), with the 
corresponding Moreau proximal mapping implemented by 10 
iterations of Chambolle's algorithm [14]. 

Table VII shows the number of calls, number of iterations. 



TABLE VII 

MRI reconstruction - Comparison 



Algorithm 


Calls to B,B" 


Iterations 


time (seconds) 


MSE 


NESTA 


1228 (1161/1261) 


307 


15.50 


9.335e-6 


C-SALSA 


366 (365/368) 


122 


12.89 


2.440e-6 



and CPU times, while Figure 4 plots the evolution of the 
objective function and constraint over time. Figure 3(c) shows 
the estimate obtained using C-SALSA (the estimate NESTA 
is, naturally, visually indistinguishable). Again, we may con- 
clude that C-SALSA is faster than NESTA, while achieving 



TABLE VI 

Image deblurring using TV - Computational load 



Expt. 


Avg. calls to B, B^ (min/max) 


Iterations 


CPU time (seconds) 




NESTA 


C-SALSA 


NESTA 


C-SALSA 


NESTA 


C-SALSA 


1 


7783 (7767/7795) 


695 (680/710) 


1945 


232 


311.98 


62.56 


2A 


7323 (7291/7351) 


559 (536/578) 


1830 


150 


279.36 


38.63 


2B 


6828 (6775/6883) 


299 (269/329) 


1707 


100 


265.35 


25.47 


3A 


6594 (6513/6661) 


176 (98/209) 


1649 


59 


250.37 


15.08 


3B 


5514 (5417/5585) 


108 (104/110) 


1379 


37 


210.94 


9.23 
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(b) 

Fig. 2. Image deblurring (uniform blur) with TV regularization: (a) evolution 
of the objective function over time; (b) evolution of the constraint ||Bx — y | 
over time. 



comparable values of mean squared error of the reconstructed 
image. 

Hi^h Dynamic Range TV Reconstruction: A related ex- 
ample that we will consider here is the reconstruction of 
images composed of random squares, from their partial Fourier 
measurements, with TV regularization (see [6], section 6.4). 
The dynamic range of the signals (the amplitude of the 
squares) varies from 20 dB to 80 dB. The size of each image is 
128 X 128, the number of radial lines in the DFT measurement 
mask is 27 (corresponding to m/n « 0.2), and the Gaussian 
noise standard deviation is tr = 0.1. 

Figure IV-C shows the original image with a dynamic 
range of 40 dB and the estimate obtained using C-SALSA. 
Figure 6 shows the evolution over time of the objective and the 
error constraint for C-SALSA and NESTA, while Table VIII 
compares the two algorithms with respect to the number of 
calls to A and A^, number of iterations, CPU time, and MSB 
obtained, over 10 random trials. It is clear from Table VIII 
that C-SALSA uses considerably fewer calls to the operators 
A and than NESTA. 



D. Image Inpainting 

Finally, we consider an image inpainting problem, as ex- 
plained in Section III-C. The original image is again the 
Cameraman, and the observation consists in losing 40% of 
its pixels, as shown in Figure 7. The observations are also 
corrupted with Gaussian noise (with an SNR of 40 dB). 




(a) 




(b) 




(c) 

Fig. 3. MRI reconstruction: (a)128 X 128 Shepp Logan phantom; (b) Mask 
with 22 radial lines; (c) image estimated using C-SALSA. 



The regularizer is again TV, implemented by 10 iterations of 
Chambolle's algorithm. 

The image estimate obtained by C-SALSA is shown in 
Figure 7, with the original also shown for comparison. The 
estimate obtained using NESTA was visually very similar 
Table IX compares the performance of the two algorithms, 
and Figure 8 shows the evolution of the objective function for 
each of them. 

V. Conclusions 

We have presented a new algorithm for solving the con- 
strained optimization formulation of regularized image re- 
construction/restoration. The approach, which can be used 
with any type of convex regularizers (wavelet-based, total 
variation), is based on a VS technique which yields an 
equivalent constrained problem. This constrained problem is 
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TABLE VIII 

Image reconstruction (high dynamic range) using TV - Computational load 



Dyn. range 


Avg. calls to B, (min/max) 


Iterations 


CPU time (seconds) 


MSE 


(dB) 


NESTA 


C-SALSA 


NESTA 


C-SALSA 


NESTA 


C-SALSA 


NESTA 


C-SALSA 


20 


1213 (1169/1273) 


226 (224/227) 


303 


76 


8.99 


7.24 


0.00241743 


0.000543426 


40 


991 (961/1017) 


227 (224/227) 


248 


76 


7.34 


7.002 


0.00432206 


0.000651107 


60 


731 (721/737) 


282 (281/284) 


183 


95 


4.92 


8.35 


0.005294 


0.00072848 


80 


617 (613/617) 


353 (350/353) 


154 


118 


4.16 


10.72 


0.00702862 


0.000664638 



1300 
1200 




2 4 6 S 10 1 


14 16 


(a) 






B £ 

C-SALSA 











(b) 

Fig. 4. MRI reconstruction with TV regularization: (a) evolution of the 
objective function over time; (b) evolution of the constraint ||Bx — y|| over 
time. 



-C-SALSA 

NESTA 



4 6 
seconds 



(a) 



-C-SALSA 
NESTA 



4 6 
seconds 



I 



(b) 

Fig. 6. TV image reconstruction (dynamic range = 40 dB): (a) evolution 
of the objective function over time; (b) evolution of the constraint ||Bx — y || 
over time. 

TABLE IX 
Image inpainting: Comparison. 





Calls to B.B-" 


Iterations 


time (seconds) 


MSE 


NESTA 


403 (401/405) 


101 


10.29 


81.316 


C-SALSA 


143 (143/143) 


47 


12.97 


75.003 



(a) 



(b) 



Fig. 5. TV image reconstruction: (a) Original image with dynamic range 
= 40 dB; (b) Estimate using C-SALSA. 



then addressed using an augmented Lagrangian method, more 
specifically, the alternating direction method of multipliers 
(ADMM). Our algorithm works for any convex regularizer 
for which the Moreau proximal mapping can be efficiently 
computed, and is therefore more general purpose than some 
of the available state-of-the-art methods which are available 
only for either (i- and/or TV regularization. Experiments on 
a set of standard image recovery problems (deconvolution, 
MRI reconstruction, inpainting) have shown that the proposed 



algorithm (termed C-SALSA, for constrained split augmented 
Lagrangian shrinkage algorithm) is usually faster than pre- 
vious state-of-the-art methods. Automating the choice of the 
value of the parameter /it remains an open question. 
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